/*
 * parmerge.c
 *
 *  Created on: Mar 12, 2019
 *      Author: sascha
 */

#include <stdio.h>
#include <stdlib.h>
#include <omp.h>

#include "merge.h"

#define MIN(x, y) (x) < (y) ? (x) : (y);
#define MAX(x, y) (x) > (y) ? (x) : (y);

void corank(int i, double A[], int n, int *j, double B[], int m, int *k) {
    *j = MIN(i, n);
    *k = i - (*j);
    int jlow = MAX(0, i - m);
    int klow = 0;
    int d = 0;
    int done = 0;
    //printf("%d %d %d %d %d %d initial\n", d, i, *j, jlow, *k, klow);

    do {
        if ((*j) > 0 && (*k) < m && A[(*j) - 1] > B[*k]) {
            d = (1 + (*j) - jlow) / 2;
            klow = *k;
            *j -= d;
            *k += d;
        } else if ((*k) > 0 && (*j) < n && B[(*k) - 1] >= A[*j]) {
            d = (1 + (*k) - klow) / 2;
            jlow = *j;
            *k -= d;
            *j += d;
        } else
            done = 1;
        //printf("%d %d %d %d %d %d\n", d, i, *j, jlow, *k, klow);
        //getc(stdin);
    } while (done == 0);
}

void merge(double a[], long n, double b[], long m, double c[]) {
    int p = omp_get_max_threads();

    int coj[p + 1];
    int cok[p + 1];

    #pragma omp parallel shared(coj, cok)
    {
        p = omp_get_num_threads();
        int i = omp_get_thread_num();

        coj[i] = cok[i] = 0;

        corank(
            i * (n + m) / p,
            a, n, &coj[i],
            b, m, &cok[i]);

        #pragma omp single
        {
            coj[p] = n;
            cok[p] = m;
        }

        int h = i * (n + m) / p;

        int j = coj[i];
        int jj = coj[i + 1];
        int k = cok[i];
        int kk = cok[i + 1];
        while (j < jj && k < kk)
            c[h++] = (a[j] <= b[k]) ? a[j++] : b[k++];

        while (j < jj)
            c[h++] = a[j++];

        while (k < kk)
            c[h++] = b[k++];
    }
}

void merge2(double a[], long n, double b[], long m, double c[]) {
    // replace this by a parallel merge algorithm
    //seq_merge1(a, n, b, m, c);

    #pragma omp parallel
    {
        int p = omp_get_num_threads();

        int j1, j2, k1, k2;
        #pragma omp for private(j1, j2, k1, k2) nowait
        for (int i = 0; i < p; i++) {
            corank(
                i * (n + m) / p,
                a, n, &j1,
                b, m, &k1);

            corank(
                (i + 1) * (n + m) / p,
                a, n, &j2,
                b, m, &k2);

            int h = i * (n + m) / p;

            while (j1 < j2 && k1 < k2)
                c[h++] = (a[j1] <= b[k1]) ? a[j1++] : b[k1++];

            while (j1 < j2)
                c[h++] = a[j1++];

            while (k1 < k2)
                c[h++] = b[k1++];
        }
    }
}

